Unquenched domain wall quarks with multi-bosons 



I. Montvay 
Deutsches Elektronen-Synchrotron DESY 
(N ; Notkestr. 85, D-22603 Hamburg, Germany 

O 

(N 

< 



Abstract 

The numerical simulation of domain wall quarks with the two-step multi- 
psj ■ boson (TSMB) algorithm is considered. The inclusion of single quark flavours, 

; as required for strange quarks, is discussed. The usage of computer memory 

can be kept relatively low, independently of the order of polynomial approxi- 
mations. Tests are performed with two flavours (Nf = 2) of degenerate quarks 

CSj . near the N t = 4 thermodynamical cross over. 

O 

1 Introduction 

Domain wall fermions JTJ, |3|, 0fl offer the possibility of improving chiral symmetry of 
lattice discretizations for fermionic theories by tuning the action parameters in an extra 
(fifth) dimension. It can be shown || that in the limit of vanishing lattice spacing in 
the fifth dimension the domain wall formulation is equivalent to the overlap formulation 
j7|, |8| which fulfills the Ginsparg- Wilson relation || for lattice chiral symmetry |L(J. The 
price of the chiral symmetry at non-zero lattice spacing is the extra dimension enlarging 
the number of degrees of freedom and, from a technical point of view, the extensions of 
the fermion matrix. It is an interesting question how much numerical simulations with 
light domain wall quarks become slower than, say, with "unimproved" Wilson fermions. 

Up to now unquenched domain wall fermions have been treated either by the Hybrid 
Monte Carlo or, in case of an odd number of fermionic flavours, by the Hybrid 



Molecular Dynamics R-algorithm [12]. (For a few examples of these simulations see, for 



instance, JT3], |TJ, [T3|, |TB[. For a recent review see |T7[ .) In the present paper the application 
of the two-step multi-boson (TSMB) algorithm |18| for domain wall fermion simulations 



is considered. This algorithm is applicable for any number of fermion flavours and has a 
tolerable slowing down towards light fermions. 

The plan of this paper is as follows: in the next Section the lattice action for two 
degenerate domain wall quarks is formulated. In Section |3] the TSM B algorithm is briefly 
recapitulated and the generalization to the case of an odd number of domain wall fermion 
flavours is discussed. In Section [| the results of test runs on 8 3 • 4 lattices near the 
thermodynamical cross over are given. 



2 Lattice action 

In this paper the domain wall fermion action is constructed according to Shamir's pre- 
scription ||. Therefore the fermion field is defined on a five dimensional hypercubic lattice 
and the light chiral fermion modes are located on two boundaries of the fifth dimension. 
The gauge field links depend only on the coordinates of the four-dimensional space-time. 
The bosonic Pauli-Villars fields subtracting the heavy fermion modes are introduced as 
in 



0- 



The complete lattice action is given by 

S = S G [U] + S F ^, % U] + S PV [&, $, U\ . (1) 
Here the standard Wilson action Sq for the SU(iV c ) gauge field U is a sum over plaquettes 

S G = /3Wl--^-ReTrfv) , (2) 

pi V iV c / 

with the bare gauge coupling given by (3 = 2N c /g 2 . In particular, in QCD the number 
of colours is N c = 3. The Pauli-Villars action in (HD Spy will be discussed later and the 
fermion action Sp for a single fermion flavour is given by 

S F = V(x',s')D F (x',s';x,s)V(x,s). (3) 

x,s; x',s' 

The four-dimensional space-time coordinates are denoted by x, x' and the fifth coordinates 
are s,s'. The domain wall fermion matrix Dp is constructed from the standard four- 
dimensional Wilson fermion matrix 

1 4 

D(x', x) = 5 X > X (4 - am ) - -^ [8 x ', x+fl (l + l^)U Xfl + 5 X > + ^ X (1 - 7/^*4'^ ■ ( 4 ) 

The notations are standard: a is the (four-dimensional) lattice spacing and fx denotes the 
unit vector in direction fx. The bare mass — m is chosen to be negative and should be 



tuned properly for producing the light boundary fermion state. In an s-block notation 
the domain wall fermion matrix is 

(a + D -aP L 

-<tP r a + D -aP L 

-aP R a + D 
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V am f P L 
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(5) 



Here rrif denotes the bare fermion mass of the light boundary fermion, P R = ~(1 + 
7 5 ), P L — ^(1 — 75 ) are chiral projectors and a = a/a s determines the lattice spacing in 
the fifth dimension a s relative to a p9| . 

The domain wall fermion matrix Dp is non-hermitean but it satisfies the relation 



D F = 75# 5 D F 7 5 i?5 



(6) 



with the reflection in the fifth dimension {R^) s ^ s = 5 s '^n b +i-s, (1 < s < N s ). This relation 
implies that the determinant of Dp is real and Dp = j^R^Dp is hermitean. Using the 
hermitean Wilson fermion matrix D = j^D the hermiticity of Dp in nicely displayed in 
an s-block form: 
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(7) 

The Pauli-Villars action is designed to cancel the contribution of the heavy fermions 
in the large iV s limit. It has a Gaussian form resulting after integration in the inverse of 
the fermion determinant. In order that the Gaussian integrals be well defined let us first 



only consider the case of two degenerate fermion flavours when the Pauli-Villars action is 



Spy = &(x", s") D F (x", s"; x', s') amf=1 D F (x', s'; x, s) amf=1 $(x, s) 



x,s: x's'\x",s" 



As it is shown by the notation, the bare mass parameter arrif is fixed here at arrif = 
1. After performing the Grassmannian integrals over the fermion fields \l/ and the 
Gaussian integrals over the Pauli-Villars fields $t, $ the result, for two degenerate fermion 
flavours, is the ratio of determinants 

det(D F )det(D F ) _ det(D 2 F ) 

(9) 



^( D F,am f =l D F,am f =l) det (L>| )0m/=1 ) 

Here we used the fact that det(-Dp) is real and det (Dp) = det(D F ). 

According to (§) the effective gauge action describing two degenerate flavours of 
fermions is a function of the squared hermitean fermion matrix Dp. The same is also 
true in case of any number of flavours, as it will be discussed in the next section. The 
non-zero matrix elements of D 2 F are, in an s-block form: 

(D%)., B = 2a 2 + 2aD r + D 2 + 5 sA (a 2 m 2 f - a 2 )P L + 5 SiNs (a 2 m 2 f - a 2 )P R , 

(D F ) SjS+ i = (D F ) s+ x jS = -a 2 - aD r , 

(D 2 F ) hNs = (D 2 F ) NsA = am f (a + D r ) , (10) 
where D r = |(7s-D + -D75) contains the Wilson term in the Wilson-fermion action: 

1 4 

D r (x',x) = 5 X > X (4 - am ) - - ^ \5 x i >x+j xU XIJt + . (11) 



/x=l 



3 TSMB algorithm for domain wall fermions 

The absolute value of the fermion determinant of an arbitrary (integer) number Nf of 
domain wall fermion flavours is, according to (Q), 

det(D F )\ Nf = {detOD*)}^ 2 • ( 12 ) 

Negative values of Nf describe the Pauli-Villars fields. (The mass parameter arrif is 
different for physical fermion flavours and for Pauli-Villars fields, but this difference does 
not play a role in what follows.) For odd numbers of flavours the sign of the determinant 
is neglected in flT2|). However, if the mass parameter is positive (rrif > 0) this sign 



is expected to be irrelevant because of the relation of domain wall fermions to overlap 
fermions which have a positive determinant if the mass is positive |5|, [(J . 

Multi-boson algorithms |2(J are based on polynomial approximations P n satisfying 

lim PJx) = x~ N f /2 (13) 

which allow to represent the fermion determinant as 

1 1 F,j det P n (D 2 ) 1 ' 

Assuming that the polynomial roots occur in complex conjugate pairs, one can write P n 



as 



P n (D 2 F ) cx H(D F - p*)(D F - Pj ) . (15) 

3=1 



This leads to the multi-boson representation 

N f n 

det(D F ) oc UdetiiDp-p^Dp-pj)]- 1 



3=1 



oc j[d<P}[d<p+} expl-i2J2<l>i4(D F - p*)(D F - Pj )] x , x cP jx \ . (16) 

[ j=l xx' ) 

Here <pj X , (j — 1,2, ... , n) are complex boson (pseudofermion) fields. 

The two-step multi-boson algorithm [18|] is based, instead of (|13j), on a polynomial 
approximation by a product of polynomials 

]\m o P£{x)P%{x)=x- N fl\ (17) 

where the first polynomial P$(x) itself is an approximation to x~ Nf ^ 2 , but it has a rela- 
tively low order. The multi-boson representation (|16|) is only used for the first polynomial 
Pi 1 ) . The correction factor P^ in 

, ~ , N f 1 

det{D F 



det Pff{D 2 F ) det P£\D 2 F ) 
is realized in a stochastic noisy Metropolis correction step with a global accept-reject 



condition, in the spirit of pl| . In order to obtain the appropriate Gaussian vector for the 
noisy correction the inverse square root of is also needed. This can be represented 
by another polynomial approximation 

P$<?)^P$V)-t ■ (19) 



A practical way to obtain P^ is to use a Newton iteration 

H% = Up!f ) + -^—] , k = 0,1,2,... . (20) 



The TSMB algorithm becomes exact only in the limit of infinitely high polynomial 
order: n 2 — * oo in (|TTD- (|T5]) . Instead of investigating the dependence of expectation values 
on ri2 by performing several simulations, it is better to fix some relatively high order n 2 
for the simulation and perform another correction in the "measurement" of expectation 
values by still finer polynomials. This is done by reweighting the configurations [^2] . The 
reweighting for general Nf is based on a polynomial approximation P^f which satisfies 

VmPW{x)P%\x)P!${x) = x~ N f/ 2 . (21) 



For more details see, for instance, J23J, [33J. 



Up to this point the particular form of the domain wall fermion action introduced in 
the previous section makes no difference compared to other applications of TSMB. The 
occurrence of negative number of flavours, as used for the Pauli-Villars fields, is the only 
new feature. However, for even number of flavours one can just use the form given by 
In this case the first polynomial for the Pauli-Villars fields is exact and no corrections 
are needed. For odd number of flavours one has to deal with polynomial approximations 
of some integer power of yfx and the machinery of polynomial approximations works as 
usual. 

Another peculiarity of domain wall fermions is that the fermion field has an extra 
index labeling the fifth coordinate. In practice this can easily lead to a situation where 
the storage of ri\ multi-boson fields in computer memory becomes a problem. (For typical 
simulation parameters including the values of ri\ see the recent studies in ^6], 27||.) 



Fortunately, in cases if the storage of the multi-boson fields is problematic, one can orga- 
nize the gauge field update in such a way that the dependence on the multi-boson fields is 
collected in a few auxiliary 3 <S> 3 matrix fields which can be easily stored in memory. The 
multi-boson fields can be kept on disk and have to be read before and written back after 
a complete boson field update. The duration of the input-output is negligible compared 
to the time of the update. 

The auxiliary 3 <S> 3 matrix fields are spin-traces over the multi-boson fields which 
can be constructed as follows. The dependence of the effective gauge field action on the 
multi-boson fields can be summarized by the formula 

S e //(£4 M ,0) = ReTr (U^SW((j>) + C^Sjgfa)) + ' " ■ (22) 
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Here U X(l is the gauge link variable to be updated. The omitted part denoted by the 
dots contain terms from the pure gauge action and other terms which do not depend on 
U xfJi . S^i 2 ^((j)) display the dependence on the multi-boson fields. gives the nearest 
neighbour contributions and S^ 2 \ which appears because of the quadratic nature of Dp, 
stands for next-to-nearest neighbour contributions. A consideration of the multi-boson 
action (|T6| ) and of the matrices in (|7|), fllCf ) gives that 



N s _ N a -1 
S xi!((p) = (fs,N s +l-s;xn + Xfs%fi) + Vffl,N s ;xn ~ ° E /s.s+l;^ 



(23) 



and 



c(2)/j.n ( f (++) r/ , - -+- f (+_) u ] +/7 f f ( ~ +) + /7 - 

^xn Vr) — Z-^ \Jx,tiu u x+ii,v i Jx,fii/ u x+[i-v,v ' u x,u Jx+i>,i/fi ' u x~u,u J x -v 
VJ^/J., u=l 

Here we used the notations 

[if = am f , 1^0 = am , x = 4 — /i + o" . 
The traces over spinor indices appearing in (JZ3"]), resp. fl2~4]) are defined as 



V/J, 
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EE Tr sp (1+7^ + 7^ + 7*7/0 
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[ 1 lvs 



i=i s=i 



(24) 



(25) 



(26) 



The indices on the multi-boson fields 4>j S ,x are as follows: j is labeling the different multi- 
boson fields as they appear in flTEp, x is the four-dimensional site and s the fifth coordinate. 
The colour and spinor indices are not shown. After performing the trace over spinor indices 
the result is, of course, a 3 <S> 3 complex matrix. 

Using the formulas (^)-(f2^) the effect of the multi-boson fields can be collected in the 
auxiliary 3 ® 3 matrices. The total number of 3 <S> 3 matrices to be stored in memory is 



4 + 3 • 12 = 40 per four- dimensional site, which is usually not a problem. The multi-boson 
fields can be updated one-by-one and kept otherwise on disk. This means that one has to 
store <pj S , x only for a single value of the multi-boson index j. Of course, due to the index 
s, the storage requirement is increasing with N s . 



4 Numerical simulation tests 

Test runs have been performed for two degenerate quark flavours Nf = 2 on 8 3 -4 lattices in 
the vicinity of the N t = 4 thermodynamic crossover. The parameter sets have been chosen 
from the points in parameter space which were investigated in Typical parameters 
were: /i = 1.9, fi f = 0.1, a = 1.0, 0.5 and 5.20 < (3 < 5.45. 

The first task is to determine the parameters of the necessary polynomials for different 
extensions of the fifth dimension N s . These will be largely influenced by the condition 
number of the squared hermitean fermion matrix Dp, because the interval [e, A] where the 
polynomial approximations are optimized has to cover the eigenvalues of Dp on a typical 
gauge configuration. A set of polynomial parameters for different N s is given in table 0. 
As it is shown by the table, the largest eigenvalues of Dp are only slightly increasing with 
N s but, in the covered range, the smallest ones decrease roughly proportional to 1/N S . 
Therefore the condition number A/e is proportional to N s . (For a discussion of bounds 
on the condition number see p8|.) 



Table 1: Parameters of the polynomials used on 8 3 ■ 4 lattice for different lattice 
spacings and extensions of the fifth dimension. The bare parameters in the lattice 
action are specified in the text. The last column shows the average number of inverter 
iterations in the quasi heatbath update step. 



o 


N s 


ni 


n 2 


n 3 


U4 


e 


A 


Iqhb 


1.0 


8 


44 


240 


300 


300 


0.011 


56.0 


5000 


1.0 


12 


56 


300 


450 


470 


0.0071 


57.0 


7900 


1.0 


16 


64 


350 


470 


500 


0.0052 


58.0 


9700 


1.0 


24 


72 


450 


640 


640 


0.0032 


59.0 


12300 


0.5 


6 


48 


270 


350 


360 


0.0071 


43.0 


6000 


0.5 


12 


64 


400 


570 


570 


0.0046 


44.0 


10100 



Table [l] refers to the case Nf = 2 which is considered in the present paper. As 

8 



discussed in the previous section, an odd number of flavours would also require polynomial 
approximations for the Pauli-Villars fields. Since the mass parameter is kept there at 
fif = 1, the corresponding condition numbers are much smaller and the polynomial orders 
substantially lower. In the runs shown by table p] the ratio of condition numbers is typically 
of the order of 10. This would require, for instance, in case of the first line of the table 
polynomial orders n\ = 32, n 2 = 60, n 3 = 90, n 4 = 100. 

A first estimate of the computation work can be given in terms of the number of 
fermion-matrix vector multiplications needed for performing a sweep over the multi-boson 
and gauge fields. An approximate formula per update cycles is: 



<IMVM 



Mqhbc^hb + Q(niN B + N G ) + 2N c (n 2 + n 3 ) . (27) 



Here Nb is the number of boson field update sweeps per cycle, Nq the number of gauge 

field update sweeps and Nq the number of noisy Metropolis accept-reject steps. It is 

assumed that a global quasi heatbath for the boson fields is performed once per cqhb 

cycles with a average number of matrix inverter iterations Iqhb- (Note that Iqhb is a 

sum over n\ inversions.) The relative contributions of the three terms in (|27|) are subject 

to optimization. The quasi heatbath is typically an important part of operations and 

Iqhb is characteristic for the total number of matrix vector multiplications Nmvm- For 

the actual values occurring in the runs see table II] which shows that, for a given value of 

i 

lattice spacing ratio a, Iqhb increases somewhat faster than Ns . 

In order to obtain an estimate for the total amount of necessary computation work, 
Nmvm has to be multiplied by the integrated autocorrelation length T int given in number 
of update cycles. The measurement of autocorrelations requires high statistics and a lot 
of computer time and is beyond the scope of the present paper. A rough estimate based 

on previous experience with TSMB tells that r int is proportional to n\. According to the 

i 

table a very rough estimate for the increase of n\ is n% oc Ni . Taking into account that a 
single matrix vector multiplication grows linearly with N s , this would finally imply that 
the increase of the necessary computation work is between cx N% and oc N%. The fast 
increase with N s favors domain wall fermion schemes with a < 1 and relatively small N s , 



as proposed in [30 



As these test runs show, the application of the TSMB algorithm for numerical simula- 
tions of domain wall quarks is possible. Simulations using alternative formulations based 
on overlap fermions, as proposed in |3l], 25] , should also be feasible along the same lines. 



The interesting question about the computation speed compared to, say, Wilson quarks 
requires a more detailed analysis including the measurement of autocorrelations. Since 
the good chiral properties of domain wall fermions develop only sufficiently close to the 



continuum limit, the performance studies have to be finally performed on large lattices. 
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